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ABSTRACT 

We investigate the dynamics of a cosmological dark matter fluid in the Schrodinger formula- 
tion, seeking to evaluate the approach as a potential tool for theorists. We find simple wave- 
mechanical solutions of the equations for the cosmological homogeneous background evo- 
lution of the dark matter field, and use them to obtain a piecewise analytic solution for the 
evolution of a compensated spherical overdensity. We analyse this solution from a 'quantum 
mechanical' viewpoint, and establish the correct boundary conditions satisfied by the wave- 
function. Using techniques from multi-particle quantum mechanics, we establish the equa- 
tions governing the evolution of multiple fluids and then solve them numerically in such a 
system. Our results establish the viability of the Schrodinger formulation as a genuine alter- 
native to standard methods in certain contexts, and a novel way to model multiple fluids. 
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1 INTRODUCTION 



On scales smaller than the Hubble distance, the evolution of struc- 
ture in the universe can be analysed using Newtonian equations for 
gravitational instability (Jeans 1928). In particular, the phase space 
evolution of the dominant cold dark matter (CDM) component is 
governed by a collisionless Boltzmann, or Vlasov, equation, cou- 
pled to a Poisson equation for gravity. When velocity dispersion 
is negligible, a fluid approximation holds, which leads to a much 
simpler representation, but one that fails to accommodate multiple 
streams in the fluid, most notably at the onset of so-called 'shell- 
crossing' in the non-linear evolution of density perturbations. 

In the linear regime, one can (by definition) describe struc- 
ture growth accurately using first-order perturbation theory (Pee- 
bles 1980) applied to the standard Newtonian equations for a self- 
gravitating fluid. In the strongly nonlinear regime, only numerical 
simulations are capable of modelling general fluid evolution and 
reproducing the structures observed in galaxy surveys; for a recent 
review, see Springel et al. 2006. However, the evolution of den- 
sity fluctuations in the so-called 'weakly nonlinear regime', those 
characteristically represented on the largest scales in galaxy cluster- 
ing surveys, have been well described by a variety of semi-analytic 
methods. These include nonlinear perturbation theory, and related 
schemes such as the Zeldovich and adhesion approximations (see 
Bernardeau et al. 2002 for a review). These methods typically ex- 
trapolate results from the linear regime into the quasi-linear regime. 

In a separate line of research, the correspondence be- 
tween the equations describing a collisionless fluid and the 
Schrodinger equation was first identified by Madelung in 1926. 
Outside cosmology, this correspondence is most commonly used in 
drawing a hydrodynamic analogy for quantum mechanics (Spiegel 
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1980; Skodje et al. 1989). Widrow & Kaiser (1993) were the first 
to apply the ideas to the problem of cosmological structure for- 
mation, using a Schrodinger field representation for collisionless 
matter, but omitting the 'quantum pressure' term to obtain a linear 
wave-equation. They developed a numerical technique to follow the 
nonlinear evolution of the field, offering an alternative to N-body 
particle mesh codes. A fuller explanation and exploration of their 
method is given by Davies & Widrow (1997). 

Coles & Spencer (2002) take a different approach, using 
the so-called Madelung transformation to establish a simpler cor- 
respondence between the self-gravitating fluid equations and a 
Schrodinger-Poisson system. This transformation means that their 
representation of the wavefunction can only be used for single 
streaming fluids, with their aim to develop a semi-analytic tech- 
nique for tracking the development of structure into the quasilinear 
regime. Short & Coles (2006a,b) develop this foundational work 
into what they term the free-particle approximation for precisely 
this purpose, in which the gravitational potential and the quantum 
pressure are neglected. They show that this approximation is use- 
ful into the mildly nonlinear regime and the resulting method is 
essentially equivalent to the well-known adhesion approximation. 
Szapudi & Kaiser (2003) develop a Schrodinger perturbation the- 
ory very similar in nature to the well-known perturbation theory of 
the standard Eulerian fluid equations. 

With the exception of Szapudi & Kaiser, all previous work in 
this area has been aimed chiefly at developing better computational 
techniques for tracking the evolution of structure, based almost in- 
cidentally on a wave mechanical formulation, which may confer 
certain advantages over rival methods; the full wave-mechanical 
system of equations that govern structure growth has not been 
a feature of interest. In this paper, our aim is to explore the 
Schrodinger formulation of the evolution of structure from a fully 
self-consistent point of view, without recourse to approximations. 
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The plan of the paper is as follows. In Section|2l we outline the 
Schrodinger formalism for fluid dynamics, summarizing the previ- 
ous work in the field and outlining our own direct approach. In Sec- 
tion[3] we obtain solutions to the Schrodinger-Poisson system for 
spatially-flat and closed homogeneous cosmological models. These 
solutions are then used in Section|4]to derive a wave-mechanical so- 
lution for a compensated spherical overdensity model. In Section[5] 
we describe an extension to the Schrodinger formalism, based on 
techniques from multi-particle quantum mechanics, that allows one 
to accommodate multiple fluids. We apply this technique to a mod- 
ified spherical overdensity model that exhibits multi-streaming by 
integrating the evolution equations numerically. Finally, our con- 
clusions are presented in Section[6] 



2 SCHRODINGER FORMALISM FOR FLUIDS 

We take as our starting point the usual Newtonian equations for 
a self-gravitating perfect fluid. The continuity, Euler and Poisson 
equations comprise the first and second momentum moments of 
the full Vlasov-Poisson system, and are exact provided the velocity 
dispersion is zero, i.e. for single streaming fluids. These equations 
are: 



- + (..V).- 
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Here p and p are the fluid velocity, density and pressure respec- 
tively, and V is the gravitational potential. In lO, we have modified 
the Poisson equation to include a term representing the repulsive 
force provided by a non-zero cosmological constant A. 

Adopting the usual description of cold dark matter as a pres- 
sureless fluid, and restricting our analysis to irrotational flow, which 
is normally considered for cosmological evolution in homogeneous 
and isotropic universe models, one can replace the Euler equation 
l|2j with the Bernoulli equation 



dt 



+ \ {V<f>f = -V, 



(4) 



with V = V<^ where (f) is the velocity potential. 

Since our equations describe a single streaming fluid, one may 
introduce the Madelung transformation 



■0 = Qe 



<.4>/v 



(5) 



In this transformation, p = — o? , so the fluid density is al- 
ways nonnegative, and p is an adjustable parameter which has the 
same dimension as 4>, namely [L^T"^]. It is then easy to show that 
the equations governing the evolution of ^ and V are 
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(6) 
(7) 



Thus, we see that the dynamics of a dark matter fluid may be 
described by the (modified) Poisson equation coupled to a nonlin- 
ear wave equation. This permits a direct formal correspondence to 
be made between the concepts of fluid dynamics and those of quan- 
tum mechanics, leading to potential conceptual advantages such as 
fluid quantities being associated with Hermitian operators. There 
are also some possible computational simplifications. In particular. 



the WKB limit for the wave-equation corresponds to the asymptotic 
study of the limit u (Spiegel 1980). 

Nonetheless, the nonlinear nature of the wave-equation {6]( 
does present some conceptual difficulties. Nonlinear Schrodinger 
equations and their applications in (non-cosmological) physics are 
treated thoroughly in Sulem & Sulem (1999) and Pang &. Feng 
(2005). The theory of nonlinear quantum mechanics is used to 
describe macroscopic quantum effects and quasi-particles such as 
solitons, which have a Hamiltonian which depends nonlinearly 
on their wavefunction. Typically used in the modelling of quan- 
tum semiconductors and superfluids, the theory differs from lin- 
ear quantum mechanics in many important ways. For example, the 
square of the wavefunction is no longer the probability of finding 
the macroscopic particle at a given point in spacetime, but instead 
gives the mass density of its constituent microscopic particles at 
that point. Operators are no longer linear, and hence wavefunctions 
no longer obey the principle of superposition; the unitary structure 
of linear quantum mechanics does not apply. 

In addition to these conceptual differences with standard quan- 
tum mechanics, the nonlinear nature of the coupled Schrodinger- 
Poisson equations (|6f and Q make them difficult to solve. The 
wave-equation has two sources of nonlinearity: the first from the 
coupling to the Poisson equation via the gravitational potential, and 
the second from the nonlinear 'quantum pressure' term (this term 
resembles a pressure gradient when interpreting quantum phenom- 
ena in terms of classical fluid behavior) 



P = 



(8) 



which acts like an effective potential. It is useful to insert a multi- 
plicative parameter fi in the nonlinear quantum pressure term in the 
wave-equation such that /3 = corresponds to leaving the term 
out, and /3 = 1 to retaining it. One then finds that the corresponding 
fluid equations remain unaltered, except for the Bernoulli equation 
which now takes the form 



f + i(v^)^ = -y + (i 



-!i)P. 



(9) 



The previous studies mentioned thus far, which apply the 
Schrodinger formalism to the dynamics of a cosmological dark 
matter fluid, typically discard the troublesome nonlinear quantum 
pressure term P. From l[8]l, we see that omitting P in the modelling 
of a dark matter fluid can be justified only if i V'l = \fp varies very 
slowly on the scales of interest. Alternatively, for fluids with pres- 
sure (such as a baryon component), one could adjust (or omit) the 
P term in the wave equation to model genuine effects of pressure 
(Coles & Spencer 2002), but this is clearly not relevant to mod- 
elling a dark matter fluid. 

We note from l[8), however, that in the limit — > 0, equation 
([9} reduces to the correct Bernoulli equation This has led sev- 
eral authors to discard the quantum pressure term P, thereby mak- 
ing the resulting wave-equation linear, and considering the result- 
ing theory only in the 'correspondence limit' v ^ ^. Although, in 
this limit, the fluid equations to which this modified Schrodinger- 
Poisson system is equivalent are those that hold physically, the limit 
^ is nonetheless meaningless in the context of the Schrodinger 
formulation itself; moreover, this assumption constrains these pre- 
vious authors to considering the study of phenomena in a certain 
limit of a parameter in their theory. Short & Coles, for exam- 
ple, tackle this difficulty by performing simulations for the low- 
est value of V possible numerically and, unsurprisingly, find their 
simulations agree best with A'^-body methods for smallest possi- 
ble V. Szapudi &l Kaiser's perturbation theory takes a more rigor- 
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ously analytic approach, but it is nonetheless based on a linearised 
Schrodinger equation, and thus is only valid in the correspondence 
limit. Furthermore, in developing their free particle approximation. 
Short & Coles rely heavily on results from linear perturbation the- 
ory, which is an entirely different mathematical representation of 
the physical process of interest. In fact, no study yet exists of the 
mathematical formulation of structure formation in the Schrodinger 
form which is not approximative in some way. 

In this paper, we therefore take a different approach and study 
the full nonlinear Schrodinger-Poisson system of equations l|6}j7]l. 
This general type of system has been studied before, but usually 
in the context of quantum mechanics, and always with A = 0. 
Spiegel (1980) was chiefly interested in discussing its correspon- 
dence to the fluid equations. Its convergence limits, well posed- 
ness, and semiclassical limits have been studied by Jungel, Mariani 
& Rial (2001), Jungel & Wang (2002) and Li and Lin (2003). Sim- 
ilar properties of simpler linear Schrodinger-Poisson systems have 
also been studied by Castella (1991) and by Abdallah, Mehats & 
Pinaud (2005). The nature of the coupling between the equations 
l|6j and ([7} means that the Schrodinger equation has what is effec- 
tively a time-dependent potential; solution methods for these have 
been studied by e.g. Park (2001) and Devoto & Pomorisac (1992). 
More recently, Erhardt and Zisowsky (2006) have explored numer- 
ical solution techniques for a spherically-symmetric Schrodinger- 
Poisson system describing the time evolution of an electron in a 
classical polar crystal. 

Adopting our more direct approach leads inevitably to more 
mathematical complication, but we believe this is more than offset 
by a clearer physical interpretation. In particular, in our approach 
the value of the parameter v is unimportant, since it does not appear 
in our analytic solutions for the wavefunction which is contrary 
to previous approaches. Also contrary to previous authors, we will 
not recast our equations in comoving coordinates, since we feel that 
this allows for a more straightforward physical interpretation of our 
results. Finally, it is important to stress that, in our cosmological 
context, the coupled Schrodinger-Poisson equations (|6}l7]( describe 
a fully classical system. Although we wish to utilise any possible 
links between quantum mechanics and the Schrodinger description 
of cosmological structure formation, caution must be used when 
drawing parallels with familiar quantum mechanical ideas, such as 
the uncertainty principle. 



of the gravitational potential V to the wavefunction which is pro- 
vided by the Poisson equation. As mentioned above, we will work 
in 'real' (i.e. non-comoving) coordinates. 

Ideally, we would like to obtain results for the most general 
case, i.e. that of a universe with arbitrary spatial curvature and a 
non-zero cosmological constant, from which all other cases can 
be derived. It is well-known, however, that such a solution would 
include elliptic functions, with which it is difficult to work an- 
alytically. Consequently, we will not consider the open universe 
case. Indeed, to maintain algebraic simplicity, we consider only a 
spatially-flat universe with a non-zero cosmological constant and a 
closed universe with zero cosmological constant. 



3.1 Spatially-flat universe 

Rather than solving the nonlinear SP system directly, in this case we 
can obtain the required solution very quickly by simply assuming 
the well-known cosmological result that the fluid density in this 
case evolves with (a globally defined) cosmic time t as 



P{t) 



Ac' ,2 / 3 



(11) 



This can then be substituted into the modified Poisson equation Q 
to obtain an equation for the gravitational potential V. The solution 
of this equation is easily found and can then be substituted into 
the (in this case linear) wave-equation ^ to obtain an equation for 
the velocity potential 0, which is also easily solved. The resulting 
wavefunction and gravitational potential are found to be 



V 



Ac2_ 
8^ 



cosech(Af) exp 



-coth(At)r 



(12) 



A „2 

— [coth^(A^)-3]r^ 



where A = ^y^A^. 

The solutions in the special case of a spatially-flat universe 
with zero cosmological constant are easily obtained from these re- 
sults by performing a series expansion of the functions p = a' and 
(f'/r'^ around A = 0, which read 



3 COSMOLOGICAL SOLUTIONS 

In this section, we discuss the solutions to the Schrodinger-Poisson 
(SP) system ||6}Q in the simple case of a homogeneous and 
isotropic cosmological dark matter fluid. In this case, the P term 
in the SP system is in fact identically zero, since — a' = p, 
as a physical observable, must also be homogeneous. Hence the 
wave-equation becomes linear without any need for approxima- 
tions. On the other hand, the velocity potential (j) and the gravita- 
tional potential V are not physical observables, and so may be non- 
homogeneous; we may only assume that they are functions of time 
t and radial distance r only. Thus, we take the ansatz a = a{t), 
(f) — 4'{t, r) in the wavefunction i/> and V = r). Consequently 
the Laplacian operator in the SP system is simply 



2d_ 

r dr 



(10) 



Although the nonlinear quantum pressure term drops out of the 
SP system, the system itself remains nonlinear due to the coupling 



1 _ c"A 
GnGf^ 24tvG 
1 c'tA 
3? 



O 



(13) 



Setting A to zero, we thus obtain the following solution to the SP 
system: 



V 



^2 



exp 



ir 



(14) 
(15) 



This result is easily checked by deriving it afresh by assuming the 
well-known result p oc in a spatially-flat universe with A = 
and following an analogous procedure to that used above. In any 
case, we note that, at at fixed t, the potential V has the form of a 
harmonic oscillator potential in r. Also of interest is the fact that 
the wavefunction itself is similar to the Green's function for the 
diffusion equation. 



© 2009 RAS, MNRAS 000.[TlfT2l 



4 Rebecca Johnston, A.N. Lasenby andM.P. Hobson 



3.2 Closed universe 



density parameters to be 



To obtain solutions for a closed universe, it is first useful to make 
a change of coordinate from cosmic time t to a 'conformal time'. 
In the usual cosmological approach, the dimensionless conformal 
time is defined as 



cAt 



(16) 



where 7? is the scale factor for the closed universe and c is the 
speed of light. In the Newtonian context in which we are working, 
however, c itself has no special meaning, and so we need to define a 
more appropriate conformal time variable. Introducing a constant, 
a, with dimensions of velocity, we can define a new dimensionless 
conformal time r) as 



T] - 



3M \ 
47r ) 



-1/3 



2/3 



dt, 



(17) 



where the constant M = [A-k /2>)pE? . For notational brevity, we 
also introduce the constant = 3M/(47r). 

In an analogous manner to our analysis of the spatially-flat 
case, rather than solving the SP system directly, we begin by as- 
suming the well-known form of the density evolution in a closed 
FRW model, as a function of conformal time rj, namely 



A' 



(18) 



where A is a constant. Following the same procedure as that used 
in the spatially-flat case, one finds that the wavefunction and grav- 
itational potential that solve the SP system are 



V = 



A 



exp 



^2^-2/3^4/3 



I afj, 

V 



-1/3^2/3 



2sin^(77/2) 



cos(j?/2) 2' 
-r 



4sin^(77/2) 

together with the following condition on the constants: 



(/iA)i/3 



SttG 
3 



(19) 
(20) 

(21) 



Indeed, one can recover the spatially-flat case (with zero cosmolog- 
ical constant) by taking the limits a ^ and A ^ in such a way 
that Ajci'' remains finite. 

It will be convenient to write the solutions ( 1191 ) and ( I20t in a 
more transparently cosmological form. To assist in this endeavour, 
it is first useful to obtain parametric equations for R and t in terms 
of the conformal time r\. Using the expressions and dlSb . and 
the relationship l l21t between the constants, one quickly finds the 
well-known results 



R 

t 



2GM .2 

— 5-sm (77/2), 

GM , . , 



(22) 
(23) 



where we have assumed the initial condition that r; = when 
t = 0. Indeed, from the first result we see that our velocity pa- 

where i?max 



is the maximum value 



rameter a = ^/2GM/Rn 
of 7? attained during the evolution. 

Using the results ( 122b and ( 123b . together with the ( 117b . dlSb 
and ( 121b . one can easily find the expressions for the Hubble and 



H 



n = 



1 dR 
'R~dt 

SnGp 



1 dr)dR _ 
R dt drj 
1 

cos2(j7/2)' 



cos(»7/2) 



2GM sin3(?7/2)' 



(24) 
(25) 



Both parameters have the correct limits at the 'big bang', rj = 0, 
and at 'turnaround', or the maximal expansion point, rj — tv. 

Using these results, it is straightforward now to express the 
wavefunction ( |19b and gravitational potential ( 120b in terms of stan- 
dard cosmological variables as 



V 



exp 



SttG 



2v 



(26) 
(27) 



It is simple to verify that these results reduce to the flat-space solu- 
tions ( 114b and J15b in the case where SI = 1 and H = 2/(3t). 



3.3 General case in terms of cosmological parameters 

It is, in fact, possible to obtain expressions for and V , in terms of 
standard cosmological parameters, that are valid for any pressure- 
less cosmological fluid. 

For 0, we have the obvious result 



0-1-1/ , 



(28) 



which follows since V0 gives the velocity, which should be Hr 
in the cosmological case. For V , we can use the Bernoulli equation 
Using the above expression for we can immediately deduce 



Since H = R/R,we have 



H + H^ = ^ = -qH\ 



(29) 



(30) 



where q — —RR/R? is the usual deceleration parameter. Thus, the 
expression for V becomes 



V 



\qH\\ (31) 



which can be seen to agree with all the individual expressions given 
above, and is also valid in the A 7^ case, for arbitrary spatial 
curvature. In that case, 

q^{\Q.^-^A) (32) 
where Q,m ~ 8itGp/{3H^) is the matter density, and Qa ~ 

In the same spirit of expressing results in terms of cosmologi- 
cal parameters, the general result for a is 



(33) 



a = constant x exp ( ^2 ' ^'^^ 



although this is less explicit than the results for V and 0. The over- 
all wavefunction in the general case is thus 

= constant x exp f ~^ / Hdt + '^^''''j ■ (3^) 
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4 MODELLING A SPHERICAL OVERDENSITY 



We now use the two cosmological solutions found in Section |3] as 
the basis for constructing an analytical wave-mechanical descrip- 
tion of the so called 'top-hat' spherical overdensity model, com- 
monly used to study nonlinear gravitational collapse (e.g. Peebles 
1980; Padmanabhan 1993). The reader should keep in mind, how- 
ever, that in all that follows we are simply describing ordinary New- 
tonian fluid evolution, albeit in a novel way. 

The model system comprises: an infinite outer fluid region of 
homogenous density; an inner fluid sphere, also of homogenous 
density, but overdense relative to the outer region; and a spheri- 
cal shell of vacuum which separates the two. fluids. The inner fluid 
sphere is usually 'compensated', that is, we can imagine forming 
the inner fluid region by removing the mass from the vacuum (gap) 
region and adding it (homogeneously) to the mass within the inner 
radius of the gap. As a consequence of Gauss' theorem (or equiv- 
alently of Birkhoff's theorem in general relativity), and because 
both the inner and outer fluid regions in this model are homoge- 
neous, each region of fluid effectively behaves as a distinct 'uni- 
verse', evolving independently of the other according to cosmolog- 
ical solutions such as those found in Section [3] The inner 'island 
universe' region remains isolated from the outer universe in which 
it is embedded via the compensation condition, which effectively 
ensures that the outer radius of the gap region evolves in such a 
way as to keep the two fluid regions distinct at all times. 

In the cosmological context, the outer region is generally taken 
to be a spatially-flat Friedmann universe, with compensation ensur- 
ing that the mean density of the system as a whole is kept at the 
critical value. The overdense inner region will evolve as a closed 
universe. It is also usual to demand that both the inner and outer 
universes have the same origin in time, i.e. that their 'big bangs' are 
synchronized. This simple spherical overdensity forms the basis of 
the 'Swiss cheese' model widely used in studies of cosmological 
structure formation. For simplicity, we will assume throughout that 
A = 0. 



4.1 Specification of the splierical overdensity model 

Our aim is to obtain a global wavefunction and gravitational po- 
tential that solve the SP equations l|6}l7]l for this physical system; 
we do this by constructing a piecewise solution from the cos- 
mological solutions derived in Section [3] and a vacuum solution. 
An illustration of the physical system is shown in Fig. [T] where 
(-i/ic:, 14), (^g, Vg) and (i/jf, Vf) represent the wave-mechanical so- 
lutions to the Schrodinger-Poisson system in the inner (closed), 
vacuum (gap) and outer (flat) regions respectively. The radii 
and Rf label the inner and outer boundaries of the vacuum region 
respectively. We will denote by V) the global wavefunction 
and gravitational potential that we construct, in a piecewise fash- 
ion, from the solutions for each region. 

Our first task is to establish a global time coordinate for the 
system. This is, in fact, straightforward, since equation ( I23t already 
links the 'ordinary' (cosmic) time coordinate t appearing in the so- 
lution for the spatially-flat outer region with the conformal time -q 
in the solution for the closed inner region. This condition estab- 
lishes the correct phasing between the two regions by setting the 
conformal time equal to zero at the zero of ordinary time. Since the 
overdense inner region is compensated precisely by the surround- 
ing vacuum region, our system evolves correctly in that the volume 
of the vacuum region vanishes as t — ^ (or r; 0). This is eas- 
ily verified as follows. The density of the flat universe varies as 




Figure 1. Illustration of the spherical overdensity. The wavefunction and 
gravitational potential that solve the SP system in the inner, vacuum, and 
outer regions, respectively, are (■0c, Vc), (i/)g, Vg), and (i/jf , Vf ). The radii 
Rc and R f label the boundaiies of the vacuum region. 



p = l/V67rGt2, and so 



1/3 



/9GM\ 



1/3 



(35) 



where M — (47r/3)pc^f is the total mass of the inner fluid region. 
This expression for R / can be compared to the equation l l22t for the 
evolution of outer radius, Rc, of the closed inner region. One finds 
that Rf/Rc is a function of conformal time rj only, and for small rj 
is given by 



Rc 20 



(36) 



0, i.e. as the big bang is 



which correctly tends to unity as rj 
approached. 

In the context of the finite inner fluid region of total mass AI 
evolving as an isolated closed universe, it is worth noting that the 
role of the velocity parameter a — ^/GM/^max, introduced in 
l ll7t . now becomes clear: it is the Newtonian escape velocity at the 
surface of the inner fluid when it reaches its maximum radius. The 
limit in which the inner fluid region becomes spatially-flat corre- 
sponds to the escape velocity a at this point tending to zero. (We 
also note that a — c would correspond to the maximum radius 
i?max being the Schwarzschild radius for a mass Af ). 



4.2 Boundary conditions 

To construct piecewise global solutions for the wavefunction tp and 
gravitational potential V describing the full system, it is necessary 
to determine the appropriate boundary conditions on these func- 
tions at the two moving surfaces Rc and Rf. 

We first consider the boundary conditions on the gravitational 
potential. At each boundary in this highly idealized system, the 
fluid 'boundary layer' is infinitessimally thin, and effectively mass- 
less. Consequently, it cannot support a force acting on it, which in 
turn implies that dV/dr, and hence also V itself, must be contin- 
uous at each boundary. The second derivative of V may, however, 
have a discontinuity, which will correspond to the gradient of the 
force being discontinuous at the boundary. This is indeed physi- 
cally correct, since in the vacuum region there will be tidal forces. 
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but clearly none in the homogeneous regions. The gravitational po- 
tential V is determined by the Poisson equation (O, with A = 0, 
which for our spherical geometry becomes 



or-' r or 



(37) 



A discontinuity in the second derivative of V thus implies that a 
should also possess a discontinuity at each boundary, as expected. 

We turn now to the boundary conditions on the wavefunction 
tj). Griffiths (1992) has previously considered the boundary con- 
ditions on a time-independent wavefunction in the presence of a 
potential possessing a discontinuity (modelled as a Heaviside func- 
tion, i.e. the derivative of a Dirac delta function). In our time de- 
pendent case, however, Griffiths' approach will not yield the appro- 
priate boundary conditions. To determine the consequences of the 
discontinuity in square-root density a at each boundary, we must 
turn to the equation governing the time development of a, which is 

da da d4> a d4> a d^cj) 
dt dr dr r dr 2 dr'^ 
For our solution V) for the full system to be consistent, the 
form for a that we know to be correct at the boundaries must satisfy 
this equation at all times. Consider, for example, the boundary r = 
Rf{t), where we make a transition from vacuum to homogeneous 
fluid. The form of a, as r increases, is 



(38) 



.{t,r) = f{t)H{r-Rf{t)) 



(39) 



where f{t) describes the time evolution of a within the outer (flat) 
region of homogeneous fluid and H is the Heaviside step function. 
Substituting this form into l l38b . one finds 

d/ fd(t) fd'^(p\ rr, 

dt r dr 2 dr^ j 



r \ or 



7?/(t)) = 0. (40) 



In order for the postulated form J39t of a to be consistent, the co- 
efficients of the step and delta functions need to vanish separately. 
This follows immediately for the coefficient of H, since on closer 
inspection this is seen to be the equation for the time development 
of a in the homogenous region, where a — f. The coefficient of 
the delta function needs to be evaluated at the moving boundary 
r = Rf{t), which means that one requires dR j /dt = d(j)/dr there. 
This is, however, exactly the statement that the boundary is moving 
at the speed which is given by the velocity potential at that point, 
and so this coefficient vanishes also. A similar analysis shows that 
this form for a is also acceptable at the other boundary. 

It now remains only to consider what, if any, conditions must 
be satisfied by the phase of the wavefunction, the velocity potential 
0, at the boundaries. The equation determining the evolution of (j) 
is simply the Bernoulli equation which in this case reads 



d<f>_ 
dt 



V. 



(41) 



Examining this equation, it becomes clear first of all that it is not 
necessary to consider the possibility of singularities in 4> at the 
boundaries: since V has continuous zeroth and first derivatives, 
then (f> will also remain singularity free. Additionally, in regions 
where ijj = i.e. the vacuum region, cj> may take any value at all, 
since it is not defined there. 

To summarize, we have shown that the only physical bound- 
ary conditions which must be satisfied for the solution ip, V for 
this system are those of continuity of the gravitational potential 



V and its first derivative. The wavefunction ij) itself need not be 
continuous at the boundaries. We have also established that our 
wave mechanical description of the system is thus far self consis- 
tent, by demonstrating that a step discontinuity in the modulus of 
the wavefunction at both boundaries satisfies the time dependent 
Schrodinger-Poisson system, as expected. 

4.3 Forming the piecewise solution 

Having established the boundary conditions on the wave mechani- 
cal description (i/i, V^) of the spherical overdensity system, we can 
now consider how we can form this wavefunction from those which 
separately describe each region. 

We first consider the wavefunction ip. The modulus of the 
wavefunction, the root-density a, is already specified everywhere, 
and thus we have no further work to do here. The phase of the wave- 
function, the velocity potential <j), is fully specified in the inner and 
outer fluid regions by ( |26| ( and J14b respectively, but it is unde- 
fined in the vacuum region since there is no fluid there. Moreover, 
our analysis above has shown that we do not require continuity of 
phase of the wavefunction at either boundary. 

We turn now to the gravitational potential V . We first note 
that, since q = in the vacuum region, the gravitational potential 
in this region must satisfy the Laplace equation, so that 



(42) 



where B and C are functions of time only and must be fixed by ap- 
propriate boundary conditions. We showed above that, for a phys- 
ically meaningful solution, we require V and dV/ dr to be contin- 
uous at each boundary. Consequently, we must match the vacuum 
solution ( I42t and its radial derivative at each boundary with the cor- 
responding solution for the inner or outer fluid region. Matching at 
the inner boundary yields the vacuum solution 



14,1 = GM 



2Rc 



whereas matching at the outer boundary gives 



g,2 



GM 



2Rf 



(43) 



(44) 



where Rc and R / are functions of time. These two results have the 
same functional form inside the vacuum, but differ by a constant. 
This is problematic: since V is defined through the vacuum, these 
solutions, unlike those for (f), do need to match. This problem is 
resolved, however, by appealing to the extra degrees of freedom 
offered to us by the fact that the velocity potential in each fluid 
region is fixed only up to the addition of an arbitrary function of 
time. From equation l l41t for the time evolution of (j), we see that for 
the vacuum solutions for V derived at each end of the boundary to 
match, we must simply adjust one or both of the velocity potentials 
at the two ends of the vacuum by functions of time alone such 
that their difference is altered by the following 



A(l> : 



3 

2Rc 



^RfJ 



(45) 



Physically, at the level of the Schrodinger equation, this corre- 
sponds to a gauge change, which alters the potential by a function 
of time, by introducing a phase term (the integral of the change in 
time) into the wavefunction. This phase term is introduced relative 
to the cosmological solution wavefunction in the relevant region. 
The solutions for individual regions can now be put together to 
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form a global solution with all boundary conditions satisfied. The 
solution J261427t is taken as the global solution (-;/), V) for r ^ .Re- 
in the vacuum region, between Rc and 7?/, the velocity potential is 
undefined, and we match V appropriately at the inner boundary so 
that the gravitational potential ( 143b holds throughout the gap. We 
then employ the gauge freedom in to adjust the velocity poten- 
tial in the outer region, r ^ i?/, so that the global gravitational 
potential is continuous. We do this by subtracting ( 145 1 from the 
phase of the wavefunction in ( I14t . which then becomes the solu- 
tion for <j) in the outer region. Equation ( 141b can then be used to 
find the corresponding gravitational potential in the outer region, 
which is correctly matched at the boundaries. The resulting global 
solution for ip - 
for r ^ Rc, Rc 
by 



= ae"*'''' and V, defined piecewise in a, (f> and V 
< r < Rf and r ^ Rf respectively, is thus given 




(46) 



undefined, 



(47) 



1 rr ^2 , 3GM 



/o 



V = 



GM 



(2I r) 



(48) 



GM 



+ 



3 



It is worth noting the physically interesting result that the 'cor- 
rection term' in the velocity potential (j) in the outer region r ^ Rf 
is proportional to the difference in the conformal time variables 
in the inner and outer regions. Also of interest is that, for similar 
physical systems, it is clear that we will not require any particular 
conditions on the wavefunction at the boundaries. As illustrated by 
our (somewhat pathological) example, step discontinuities ensure 
that the wavefunction and its first derivative certainly are not re- 
quired to be continuous there, although we can imagine situations 
where they may be; an example of this might be a smoothly varying 
density of fluid with a boundary represented by an abrupt change in 
fluid velocity at a given radius. This behaviour is in marked contrast 
to familiar systems in quantum mechanics, where it is the gravita- 
tional (or other) potential which generally is considered to have dis- 
continuities, but the wavefunction must satisfy physical boundary 
conditions. 



4.4 Comparison with quantum Schrodinger systems 

It is interesting to consider in more detail how our results compare 
with what may be considered a similar system in standard quan- 
tum mechanics. Surprisingly, the problem of which boundary con- 
ditions to apply at a moving boundary is a matter of debate in the 
literature: an example of the confusion is provided by Samura and 
Ohmukai (2006). The authors consider reflection and transmition 
of an incident wavefunction at a moving potential step using the 
standard time-dependent Schrodinger equation. 

They claim that the boundary conditions one would expect to 
apply, namely that the wavefunction and its first derivative are con- 



V(x) 



Vo 



Vincident 
► 


Vtrarsmitted 
^ 


< 

^reflected 











Figure 2. Reflection and transmission of a wavefunction at a boundary of 
height Vo, moving at speed v, after Samura and Ohmukai (2006). 



tinuous, are not in fact valid. Since the position of the potential step 
is a function of t, they claim that one cannot use the continuity of 
the first derivative because in taking this condition, the boundary 
is fixed in time. Instead, they argue, one should use an alternative 
boundary condition, that the phase of the wavefunction is continu- 
ous at the boundary, and they proceed to use this to derive expres- 
sions for the reflected and transmitted waves. 

Although this condition yields the correct final results, it is 
straightforward to show that the authors are incorrect in claiming 
that the continuity of the first derivative is not a valid boundary con- 
dition to apply. One can in fact use the standard boundary condi- 
tions to derive the correct results, simply applying a Galilean trans- 
formation to the stationary barrier case. 

In this approach, we take the usual Schrodinger equations left 
and right of the boundary, but take the three wavef unctions, tljijipr, 
and i/'t as functions of x' — x—vt. We examine the case in which, if 
the step were stationary, il^t would be purely evanescent, so that we 
can take the limit to an infinite step, which is the case most closely 
resembling our cosmological Schrodinger system. By applying the 
usual boundary conditions of continuity of the wavefunction and 
its first derivative at the boundary, x = vt, one finds the following 
results for the wavenumbers fcr and kt in terms of fci (setting h = 
1): 

kr = fci - 2mv, (49) 
kt = -imv + \J -m?v'^ + 2vmki - k-? + 2Vom. (50) 

These yield wavefunctions which satisfy the Schrodinger equations 
in each region, and they reduce to the correct results in the case 
where n = 0. A series expansion of the first derivative of the wave- 
function in the case that Ya ^ oq shows that there is a finite slope 
at the boundary, demonstrating that this method works for both the 
finite and infinite potential step. 

This ambiguity in the published literature highlights the 
fact that determining which boundary conditions to apply to 
Schrodinger wavefunctions, even in the usual quantum mechani- 
cal context, is not a trivial matter. It is interesting that in the case 
of our classical gravitational system, neither the matter wavefunc- 
tion nor its derivative need be continuous, in contrast to the usual 
quantum case. Instead, we require that conditions be placed on the 
gravitational potential which is coupled to the wavefunction. 



5 MULTIPLE FLUIDS 

The fluid description of cosmological dark matter breaks down 
when multistreaming occurs, most notably at shell-crossing. In 
general, the onset of shell-crossing occurs after the onset of non- 
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linearity; siiell-crossing tiius represents a natural 'barrier' beyond 
which fluid-based methods cannot follow the evolution of structure 
into the fully nonlinear regime. Attempts to describe structure in the 
quasi-linear regime using a fluid description therefore resort to ap- 
proximate methods to extend the approach beyond shell-crossing, 
e.g. the Zeldovich, adhesion, and free-particle approximation meth- 
ods. 

The Schrodinger formalism is also based on a fluid descrip- 
tion, and consequently shares these limitations. Nonetheless, as in 
other fluid-based approaches, one can model more general physical 
systems by considering them to be made up of multiple fluids. One 
still cannot accommodate multi-streaming within each individual 
fluid, but one can model some multi-streaming situations by con- 
sidering each stream as a separate fluid. In this section, we there- 
fore extend our previous description by moving from a single fluid 
to the multiple fluid case, inspired by the application of techniques 
from multiparticle quantum mechanics. We do this by assigning a 
distinct wavefunction to distinct regions of the CDM fluid, each de- 
scribed by its own wave-equation, which interact collectively via a 
common coupled Poisson equation. 

The toy cosmological system illustrated in Fig. [T] can always 
be successfully described using the Schrodinger formulation de- 
scribed in the previous section because the inner and outer regions 
of fluid never overlap, i.e. its inner and outer 'shells' never cross. 
There is thus a single unique fluid velocity for each point in the sys- 
tem, throughout its evolution. This formulation would break down 
at multistreaming because the Madelung transformation ansatz l[5) 
for the form of the fluid wavefunction assumes a single-valued ve- 
locity field. 

A good toy model for illustrating multi-streaming in cosmol- 
ogy can be constructed, however, by considering a physical system 
similar to that discussed in Section [4] but with appropriate initial 
and boundary conditions to ensure that the two fluid regions do 
overlap at some stage. We now describe how to construct a new, 
two-fluid description for such a system that accommodates this 
multi-streaming. 



form 



ip = i/)(ri,r2,t), 



(51) 



where r\ and r2 are spatial positions for an element within fluid 1 
or fluid 2 respectively, and t is the overall common time variable; 
we have thus introduced a (6 + 1) -dimensional parameter space. 
As for the single fluid case, we spilt ^ into a magnitude and phase 
written as 

V'(ri,r2,t) = a{rx,r2,t)e^p[i(f){rx,r2,t)/v]. (52) 

We also introduce the gradient operators Vi (i = 1, 2) for each 
fluid, where we adopt the convention that bold quantities 'live' in 
the individual fluid spaces, and non-bold quantities refer to the to- 
tal space (with 6 'spatial' dimensions). Using Cartesian coordinates 
Ti — {xi, yi,Zi) in each fluid space, for example, the correspond- 
ing Laplacian operators (for i = 1,2) are 

V^ = ^ + ^ + ^ (53) 

oxf oyf ozf 

We now assert that the overall two-fluid wave-equation is 



dip 



(54) 



dt " ~T 

where Vi,2 is the common gravitational potential and the quantum 
pressure term has the form 



(v? + vi)i^l 



(55) 



This expression for P is certainly the natural generalisation to the 
two-fluid case of the equivalent P-term for a single fluid. It also 
agrees with the 'many-body quantum potential' Q given in equa- 
tion (7.1 .4) of Holland (1989), although our quantum pressure term 
is the negative of this, since here we are seeking to remove the term 
that would otherwise appear in the fluid equations of motion. 

A useful quantity relating fluid density and velocity in de 
Broglie-Bohm theory is the probability current. The natural gen- 
eralisation of the probability current to two fluids, is 



j = — [(VV')V*-^(V^)* 

Zmt 



(56) 



5.1 Multiparticle wavefunctions and multiple fluids 

Since we have previously adopted a wavefunction approach to 
modelling a single fluid, it is natural to consider whether multi- 
particle methods in quantum mechanics can assist in providing a 
wavefunction approach to multiple fluids. For simplicity, we will 
confine our attention to just two fluids. 

Several questions immediately arise, such as what is the cor- 
rect form of the quantum pressure and potential terms that should 
be used in the 'two-particle' case, and whether 'entangled' states 
have any role to play in the classical evolution of two coupled flu- 
ids, or whether factorisable states are the only ones relevant. 

In the single-fluid case, the approach we have been using is 
strongly related to the de Broglie-Bohm formulation of quantum 
mechanics (for a full account see, for example, Holland 1989). It 
is natural, therefore, to turn first to the many-particle version of 
de Broglie-Bohm theory. We now wish to adapt this approach to 
obtain an equation that describes the evolution of two fluids, i.e. 
describes their velocities and particle densities in time and space, 
using a Schrodinger formalism. Of course, we also need to show 
how one can obtain individual densities and velocities for each of 
the fluids from the resulting equation. 

In the two-fluid case, the overall wavefunction now has the 



where V = Vi + V2. The interpretation of J as a current in six 
dimensions, corresponding to a flow of 'probability' into and out of 
regions, just as for the single-fluid case, follows from 



dt 



+ V • J = 0, 



(57) 



which is easy to prove by using the two-fluid wave-equation i54l . 
Using the split ( 152b into amplitude and phase terms, we find 

J = q'(Vi<?!> + V2<^), (58) 

and then identifying the individual space velocities 

vi = Vi0, V2 = V2(^ (59) 

we have 

J = C? (Vx + V2) . (60) 
The conservation equation l |57| ( can now be written 

^ + Vr(Q'?;i) + V2-(a'«2) =0. (61) 

The other equation we require is that relating the com- 
mon gravitational potential Vi,2 to the densities of the two fluids 
sources. The term Vi,2 (for which we effectively include a copy for 
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each fluid in the wave-equation) is taken to satisfy the (single-fluid) 
equation 



VVl,2 ='17VG (PI+P2) 



(62) 



i.e. the potential under which both fluids move is generated by the 
sum of their densities. 

The wave-equation J54t and the Poisson equation i62\ give 
what seems a plausible 'two-particle' quantum description of two 
interacting classical fluids. The most basic test that any many-body 
equation must satisfy is that it works for a factored state, for which 

tp{ri,r2,t) = ?/)i(ri,t)^2(r'2,t)- (63) 

Inserting this form of into the wave-equation i54\ . dividing 
through by and applying the usual separation of variables 

argument (here to the coordinates ri and r2), we obtain the two 
individual fluid equations 



"1 1^ 2 

IV—- = -—Vxipl + Vl,2^1 

ot 2 

di>2 -rr-i I , T/ ; 

'-HT = V2l/'2 + Vl,2V'2 ■ 

ot 2 



aiv-il 
^'vil^2| 

2|^2| 



■l/'2- 



(64) 



The wave-equation l l54b thus makes sense for factored states. How- 
ever, generalizing to 'entangled states' appears problematic. This 
can be seen immediately in equation J62t for the gravitational po- 
tential, which requires knowledge of (the sum of) the separate fluid 
densities pi and p2. Given a general two-fluid wavefunction ip, it 
is not clear how these would be obtained. Given i/;, we can suc- 
cessfully form the two separate particle velocities, via the defini- 
tions ( |59t , which only require that the total phase, (j), be available. 
However, knowledge of the overall magnitude, o? — ji/ip, is not 
enough to deduce the separate densities. Thus, in the general entan- 
gled case, we have no way to generate pi and p2 separately from 
ip, and furthermore, it appears that we cannot find their sum either 
(to put into the right hand side of (I62t), so that the equation set 
becomes incomplete as it stands. 

A possible solution to this problem is to assume continuity 
in the two spaces separately, which seems a reasonable physical 
requirement. Supposing that 

t + - t 

and given initial density distributions pi and p2, there is enough 
information available to propagate the two densities individually 
forward in time, since i;i and V2 are separately available. Thus 
there is in principle an answer available, although it would not be a 
useful one numerically. 

We note further that the intuitive identification o? = pip2 
is not even necessarily valid in the non-factored case, since the 6- 
dimensional continuity equation ( |60t . plus the assumption of con- 
tinuity for each fluid separately, can be used to show that 

„2 



V-(p2W2)=0, (65) 



dt 



In- 



plP2 



- (wi-Vi lnp2 + W2-V2 In pi) 



(66) 



where the d/dt represents a comoving or streamline derivative in 
the 7-dimensional space: 



d_ 

dt 



9 

= — +t;i-Vi 

ot 



- V2- 



(67) 



In the factored case, the right-hand side of l |66t will vanish, 
meaning that the identification of with pip2 is maintained along 
the (joint-fluid) motion. However, in the non-factored case, there is 
no guarantee that right-hand side of l |66t will indeed vanish, and we 
have a further problem in interpreting the wavefunction. 



To summarize our answers to the questions posed at the begin- 
ning of this section, we now have a good two-fluid wave-equation 
description of the motion of two coupled fluids, for which factored 
states make sense and return sensible equations in each particle 
space (with the coupling still present due to the potential, but not 
explicitly at the coordinate level). The notion of non-factorisable 
states looks problematic, but given that we are attempting to de- 
scribe classical, rather than quantum, fluids, this is perhaps unsur- 
prising. 



5.2 Numerical two-fluid evolution with spherical symmetry 

We now apply the two-fluid wave-equation formalism developed 
above to model explicitly the evolution of a spherically symmet- 
ric system consisting of two physically identical, mutually self- 
gravitating fluids, each of which is described by individual wave- 
functions tpi and tp2 respectively. As previously, the wavefunctions 
are given by 



1p2 



aie 



026 



(68) 
(69) 



where ai(r, t) and <j)i{r,t) and Q2(r, t) and (j>2(r,t) are the 
square-root densities and velocity potentials for each fluid. The 
equations that these wavefunctions satisfy are given by J64t , which 
must, in general, be solved numerically. 

The simplest way to solve the wave-equations numerically is 
to take their real and imaginary parts, and to solve for the amplitude 
and phase of each wavefunction. We therefore solve the following 
system of five coupled equations for the evolution of the two fluids: 



9qi _ da\ dcj>i ai d^4>i 
dt dr dr 2 dr'^ 

2 



dt 2 I ar 



d4>i 



da2 



dcX2 d4>2 Ol2 d (j}2 

dr dr 



2 dr^ 

2 



Ql 0<pl 

r dr 



a.2 0<p2 
r dr 



(70) 



(71) 



(72) 



(73) 



(74) 



dt ^ \dr 

-TTY H TT = 47rG (ai + Q2) 

ar^ r dr 

It is worth noting that equation l |74t for the gravitational potential 
V does not contain a time derivative and hence, for fixed t, it is, in 
fact, an ordinary differential equation in r. Moreover, we note that 
the parameter v does not appear in the above system of equations. 

To solve the system of equations l |70H74t numerically for the 
five variables ai{r,t), (f)i(r,t) (i = 1,2) and V{r,i), one must 
first specify the appropriate initial and boundary conditions on the 
solution. The initial conditions at some time t = Q (say) are as- 
signed by specifying (r , 0) and (f)i{r,G) (i — 1 , 2) . At any given 
t (including t — 0), knowledge of the ai is sufficient to determine 
the gravitational potential V by solving l |74t . In this spherically- 
symmetric case, the general solution of i74\ . at a given time, is 



V = 



Ci 



-dr 



C2, 



(75) 



where ci and C2 are arbitrary functions of time. Clearly C2 is merely 
an overall additive term which fixes the zero of gravitational poten- 
tial, and ci expresses whether there is a point mass at the origin. We 
choose ci = C2 = to fix F in a simple manner. Furthermore, we 
will assume that the total fluid density is non-zero at the origin. This 
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means that the dominant behavior of the inner integrand is oc r^, 
and thus that the dominant behavior of the outer integrand is oc r. 
With these assumptions, one can perform the double integration in 
( 175 1 using a set of two interleaved Simpson's rules, which we have 
found to work to high accuracy. 

It remains to consider the boundary conditions on at and <j)i 
at each end of the spatial domain, i.e. at r = and r = rmax. We 
first focus on the conditions that must apply at the origin. There is 
a physical requirement at r = that, in the absence of a singularity 
there, it should be neither a source nor a sink. Since, as mentioned 
above, we are assuming the total fluid density ai + ai > at 
the origin, we thus require dcj>\/dr = and d(f)2/dr — (i.e. zero 
fluid velocities) there. There is also an arbitrary additive constant in 
each <j>i to fix and we do this by requiring that </>i = at the origin. 
These two boundary conditions for each <j!>i are compatible with 
the equations of motion, since dcpi jdr = and = at r = 
combine to mean that d<i)i jdt is also zero there, so </>i = will be 
maintained. A further detail that has to be addressed at r = is 
how to evolve each forward in time given that 1/r appears in 
one of the terms of da.il dt. Since <j!>i and d(j)i/dr are both zero at 
r = we can assume that the dominant behavior of (f>i is oc r^. In 
this case, we can combine the last two terms of dai /dt to obtain 
for each fluid 



dai 
'dt 



(76) 



If (j)i is in fact proportional to a higher power of r, then this formula 
is still correct, since it then returns zero. It is worth noting in the 
above prescriptions, there is no requirement that the ai have zero 
gradient at r = 0. Thus a cusp in the ai distributions at the origin 
is physically allowed, and does not mean that there is a singularity 
there. 

We now consider the boundary conditions we must apply at 
the outer edge r = rmax of the spatial domain. It is, in fact, pos- 
sible to evolve the equations ( l70M74t numerically without fixing 
boundary conditions for ai and (j)i at the outer end point of the spa- 
tial grid, provided we take r = rmax to be sufficiently large that the 
total fluid density is always zero there. However, any increase in the 
size of the spatial domain comes at a heavy computational cost for 
a fixed spatial resolution, so we choose to apply suitable boundary 
conditions at the outer end point of the spatial grid. To do so, we 
use the 4 rightmost end points of the grid to calculate third-order 
spatial derivatives of ai and (/>i at r = rmax and thereby predict to 
third-order in r the values of these variables at an 'extra' endpoint, 
one grid step beyond the desired boundary. Since such a scheme is 
accurate to third order, it is sufficient for the consistent solution of 
our system of second-order PDEs. 

Having specified the above initial and boundary conditions on 
the solution, the equations | |70I474| ) can then be integrated numeri- 
cally. The most straightforward numerical approach is to use a sim- 
ple first-order time-stepping technique in which the initial distribu- 
tions ai (r, 0) and (j)i{r,0) are 'marched forward' : at each time-step 
t 1-^ t + At, we calculate the required spatial partial derivatives 
on the right-hand sides of ( |70ti73t by taking differences and then 
add At X dai/dt and At x d(j)i/dt to ai and 4>i (respectively) 
at each grid point. It is well-known, however, that such a scheme 
is prone to numerical instability (a fact that we unfortunately ver- 
ified when we implemented it). We therefore instead employed a 
fourth-order Runge-Kutta integration technique (see e.g. Press et 
al. 2007), which calculates spatial derivatives (again by simple dif- 
ferencing) at several points throughout the time interval through 
which the solution is being advanced. This offers a significant in- 



crease in accuracy over a standard Euler, or simple finite differ- 
encing, technique, which calculates spatial derivatives only at the 
beginning of the time interval. 



5.3 Toy physical system 

To demonstrate that our approach can indeed be used to describe 
the evolution of a (spherically-symmetric) two-fluid system after 
the fluids overlap, we apply it to a simple physical system similar 
to the spherical overdensity model considered analytically in Sec- 
tion |4] In that case, however, the inner and outer fluid regions did 
not overlap, the outer fluid region was of infinite spatial extent, and 
the density profiles possessed sharp step-function discontinuties. 
We consider here instead a modified physical system, as follows. 
First, we consider both inner and outer fluid regions to be of finite 
spatial extent. Second, we specify initial density and velocity distri- 
butions for both fluids that are consistent with flat-universe evolu- 
tion, since it is straightforward to construct such a system in which 
overlap of fluids is guaranteed to occur Finally, since our numerical 
solution scheme relies on obtaining accurate spatial derivatives, it 
is clearly feasible only for initial density profiles that are smoothed 
to some degree. Nonetheless, the smoothing is sufficiently slight 
that we can compare the results we obtain numerically to known 
analytic results for the evolution of an annulus of a flat universe, as 
found in Section|4] 

The assumed initial ai and <j)i distributions for each fluid are 
shown in Fig. [3] together with the resulting initial gravitational po- 
tential V derived from M5\ using a single application of the Simp- 
son's rule method described above. The initial ai profiles are step 
functions smoothed with a relatively narrow tanh function. The 
initial profile for the inner fluid is that of a flat universe, i.e. 
quadratic in r, and the velocity potential for the outer fluid was 
taken to be 4>2{r, 0) — 0.01 — O.lr^. We set the outer boundary at 
rmax = 5, and used 401 grid points to sample the spatial domain. 
The equations were evolved over 651 time steps, with a total inte- 
gration time of 3.5. It was checked that the ratio of spatial and time 
resolutions Aa; / At remained greater than or equal to the greatest 
fluid velocity in the system throughout the entire integration time. 

It is worth mentioning how the accuracy of the numerical in- 
tegration is affected by the choice of the initial conditions. As ex- 
pected, very sharp edges on the ai profiles lead to inaccuracies 
in the calculations spatial derivatives via finite differencing when 
the 'width' of the profile edges is close to the spatial grid-spacing; 
these inaccuracies can then propagate into the rest of the spatial do- 
main. Also, if either fluid velocity is too fast at any point, such that 
it exceeds Ax/ At, then inaccuracies again propagate through the 
numerical solutions. This may occur directly by setting the initial 
velocity too high, but high velocities can also occur later in the evo- 
lution if the initial density of either fluid is set very large, resulting 
in faster gravitational collapse. 

Given the initial conditions illustrated in Fig. [3] the resulting 
ai and (fii distributions for each fluid after evolving the system 
to t = 3.5, together with resulting gravitational potential V, are 
shown in Fig. |4] We see that, as the two fluids overlap, new fea- 
tures develop in the profiles. In particular, we note a pronounced 
caustic in fluid 2 at a finite distance from the origin, which signals 
that fluid is 'piling up' at this radius as a result of shell-crossing. 
Fluid 1 also possess a corresponding 'hump' at radii larger than the 
position of the caustic in fluid 2. 

Finally, to check our code, we also considered a model in 
which the initial density distribution of fluid 1 was set to zero ev- 
erywhere, but the initial density of fluid 2 has the same profile as in 
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Figure 3. The initial square-root densities, a\ and «2 (left), initial velocity potentials, 4>\ and 4>2 (middle) and the initial gravitational potential V (right) for 
the numerical integration of the spherically-symmetric two-fluid system. The solid lines refer to the fluid 1 and the dashed lines to fluid 2. The outer fluid (fluid 
2) is given a small inwards initial velocity, as shown. 




Figure 4. As for Fig. [3] but after evolving the system to a time t = 3.5. 



Fig. [5] (left panel) and its velocity potential was taken to be that of 
a spatially-flat universe. In this case, aside from minor differences 
resulting from the smooth edges of the density distribution of fluid 
2, we expect this fluid to evolve as (an annulus of) a spatially-flat 
universe (with A = 0), for which the analytic solution is given by 
il4i and il5l . Our numerical results were found to agree well with 
this analytic solution. This example also demonstrates that, because 
of the (slight) smoothing of the initial density profile of fluid 2, the 
total density at the origin r = was non-zero to machine precision, 
thereby satisfying our earlier assumption in solving equation ( 175 1 
for the gravitational potential. Hence our scheme can still be used 
in many cases to evolve systems which 'appear' to have a hole at 
the centre. 



6 CONCLUSIONS 

We have considered the use of the Schrodinger formalism for fluid 
dynamics to model the evolution of cosmological dark matter. Our 
approach differs from previous work in this area in that we consider 
the full non-linear Schrodinger-Poisson (SP) system, without mak- 
ing any approximations. In particular, we retain the 'quantum pres- 
sure' term, which previous authors have discarded to regain a linear 
wave-equation. By so doing previous authors have had to work in 
the so-called correspondence limit u ^ 0, where ;^ is a parameter 
that changes the spatial and velocity resolution of the solutions. 

We have shown that cosmological solutions to the full SP sys- 
tem can be readily obtained for closed and spatially-flat universes. 
Moreover, these solutions can be combined to obtain a global solu- 
tion to the SP system for a spherical overdensity. This necessitated 
determining the appropriate boundary conditions on the solutions 
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for each homogeneous region of the model, which differ from those 
usually assumed in quantum mechanics. We also highlight that the 
issue of what conditions to apply to wavefunctions at a moving 
boundary remains a matter of debate in the quantum mechanics lit- 
erature. 

We consider how the Schrodinger formalism for classical flu- 
ids can be extended to multiple fluids and derive an appropriate 
wave-equation for the two-fluid case, which could be easily gener- 
alised to more fluids. In particular, we find that for consistency of 
the equations the multi-fluid wavefunction must be factorisable into 
wavefunctions for each fluid separately; 'entangled' states appear 
to be problematic. Nonetheless, our approach provides a new and 
consistent method for dealing with multiple fluids. We illustrate our 
method for describing multi-streaming fluids by applying it to the 
evolution of a modified spherical overdensity toy model, in which 
the two fluids eventually overlap. We integrate the SP equations nu- 
merically to obtain physically sensible results, which illustrate the 
formation of cusps. 

We consider our investigations to show that the full non-linear 
SP system has the potential to be a useful tool for theoretical in- 
vestigations of cosmological fluid evolution. The non-linearity of 
the equations can be troublesome for analytical investigations (al- 
though not always), but poses few problems for following the fluid 
evolution numerically. In any case, we believe these difficulties are 
offset by the advantages it offers in ensuring that the fluid density 
is always positive. In future work, we plan to investigate linear per- 
turbation methods within the Schrodinger formalism to determine 
whether such an approach offers any advantages over traditional 
first-order perturbation methods applied to the standard fluid equa- 
tions in following cosmological structure formation in the linear 
regime. 
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